Innovalab
Introduccion a INLA

1 Analisis Predictivo

El análisis predictivo emplea datos históricos para predecir eventos futuros. Normalmente, los datos históricos se utilizan para crear un modelo que capture diversos patrones temporales. Este modelo predictivo se usa entonces con los datos actuales para predecir lo que acontecera. Existen diversas metodologias para llevar acabo predicciones, en este taller optamos por el enfoque bayesiano para llevar acabo la prediccion.

1.1 Enfoque bayesiano : Distribucion predictiva aposteriori

El enfoque bayesiano, en el contexto del analisis predictivo, enfoca su analisis en el calculo de la distribucion predictiva aposteriori. Esta distribucion emplea la informacion de los datos asi como la distribucion aposteriori de los parametros de interes . Empleando estos insumos la distribucion predictiva podemos generar datos futuros; es decir, datos faltantes en un determinando horizonte.

\[\begin{align*} f(y_{pred}|y_{obs}) = \int_{\theta} f(y_{pred},\theta|y_{obs})d\theta =\int_{\theta} f(y_{pred}| \theta y_{obs})f(\theta|y_{obs})d\theta \end{align*}\]

Actualmente existen diversas metodologias para llevar acabo el calculo numerico de esta distribucion las cuales en R corresponden al empleo de un determinado paquete en. Algunos paquetes ampliamente conocidos son los siguientes :

  • STAN
  • Winbugs
  • Jags

Relativamente hace no mucho tiempo la metolodogia de la Aproximación anidada integrada de Laplace o INLA(por sus siglas en ingles) ha cobrado importancia por su relativa eficiencia para estimar modelos complejos.

1.2 Prediccion como datos faltantes

El analisis de prediccion como en la practica puede ser entendida como predecir datos desconocidos en un determinado horizonte. En ese sentido estos datos a futuro pueden ser interpretados como datos faltantes o missings values, los cuales son el objeto de la predicción.

Año Semana Numero de muertes
2019 1 0
2019 2 1
2019 3 4
2020 1 0
2020 2 1
2020 3 4
2021 1 NA
2021 2 NA
2021 3 NA

Para llevar acabo predicciones empleando INLA, una aproximacion es justamente la anterior. Una vez dispuestos los datos como en la tabla anterior, basta con ajustar un modelo tipico en INLA, ya que INLA internamente calculara la distribucion predictiva aposteriori .

2 Analisis Descriptivo

library(tidyverse)
library(INLA)
library(yardstick)
library(gt)
library(innovar)
library(spdep)
library(reshape2)

db       <- readRDS("db_excess_proc.rds")
db.lima <- db %>% 
            filter(prov=="LIMA")


db.lima %>% 
ggplot()+
geom_line(aes(x=date,y=n),color="#011f4b")+
geom_point(aes(x=date,y=n),color="black")+
facet_wrap(vars(sexo))+
ylab("Numero de muertes")+
xlab("")+
ggtitle("Evolucion del Numero de muertes por género")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5,size = "22"),
      strip.background =element_rect(fill="#011f4b"),
      strip.text = element_text(color = "white",face = "bold",size = "5pts")) 

2.1 Modelamiento con INLA

Para llevar acabo el ajuste de modelos empleamos la funcion inla(), la cual tiene diversos parametros, algunos de los mas importantes son :

  • data : Un objeto tipicamente de la clase data.frame el cual representa los datos para ajustar el modelo
  • formula : Un objeto de la clase formula que especifica la ecuacion que pretendemos ajustar como por ejemplo inla(target ~ 1 + x1). En la formula podemos especificar efectos lineales (introduciendo el nombre de la variable) o no lineales (empleando f()).
  • verbose : Una variable del tipo boolean, la cual señala si se desea mostrar en consola el proceso de convergencia
db.femenino   <- db                                   %>%  
                 filter(sexo=="FEMENINO")             %>% 
                 mutate(n=ifelse(annus == 2019,NA,n))



db.masculino  <- db                                   %>%  
                 filter(sexo=="MASCULINO")            %>% 
                 mutate(n=ifelse(annus == 2019,NA,n))      


db.n.femenino   <- db.lima                            %>%  
                   filter(sexo=="FEMENINO")           %>% 
                   mutate(n=ifelse(annus == 2019,NA,n))  


db.n.masculino  <- db.lima                             %>%  
                   filter(sexo=="MASCULINO")           %>% 
                   mutate(n=ifelse(annus == 2019,NA,n))
lima.m.m0 <- inla(n ~ 1 + annus,
                  verbose         = F,
                  data            =  db.n.masculino
                  )

lima.f.m0 <- inla(n ~ 1 + annus,
                  verbose         = F,
                  data            =  db.n.femenino
                  )

Los parámetros antes detallados son los esenciales para ejecutar el ajuste de un modelo empleando INLA.Sin embargo, algunos parametros extra a tener en consideracion son los siguientes :

  • family:Objeto de clase character.Este parametro es crucial, pues determina la distribucion de la variable objetivo, por defecto se encuentra en family=Gaussian.
  • control.compute:Objeto de la clase list permite especificar el calculo de criterios de informacion tales como aic,dic,waic.
lima.m.m0 <- inla(n ~ 1 + annus,
                  verbose           = F,
                  data              =  db.n.masculino,
                  control.compute   = list(dic=TRUE, waic=TRUE,cpo=TRUE),
                  control.predictor = list(link = 1),
                  )

lima.f.m0 <- inla(n ~ 1 + annus,
                  verbose         = F,
                  data            =  db.n.femenino,
                                    control.compute   = list(dic=TRUE, waic=TRUE,cpo=TRUE),
                  control.predictor = list(link = 1),

                  )

2.2 Prediccion: Enfoque Temporal

En el analisis de prediccion (forecast) tipicamente se emplea datos en el formato de series de tiempo para llevarla acabo. Para ello comunmente se asume que las series de tiempo pueden ser descompuestas del siguiente modo :

\[\begin{align*} y_{t} = S_{t}+T_{t}+e_{t} \end{align*}\]

donde :
\(S_{t}\) : Estacionalidad
\(T_{t}\) : Tendencia
\(e_{t}\) : Ruido blanco

2.2.1 Ajuste y prediccion

Los componentes antes detallados pueden ser modelados en el contexto de INLA empleando distribuciones apriori para cada componente tales como los siguientes:

  • linear : variable
  • AR1 (Proceso autoregresivo de 1er orden) : f(variable, model = "ar1")
  • RW1 (Paseo aleatorio de 1er orden) : f(variable, model = "rw1")
  • RW2 (Paseo aleatorio de 2do orden) : f(variable, model = "rw2")

Dado el comportamiento que exibe la series que prentedemos modelar, para los componentes :

  • Tendencia : Emplearemos una distribucion apriori del tipo autoregresivo(AR1) o linear
  • Estacionalidad : Emplearemos una distribucion apriori de tipo de paseo aleatorio de 1er o 2do orden

2.2.1.1 Hombres

#inla.list.models("likelihood")

lima.m.m1 <- inla(n ~ 1 + f(annus,model = "ar1"),
                  control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
                  control.predictor = list(link = 1),
                  verbose         = F,
                  family          = "nbinomial",
                  data            =  db.n.masculino
                  )

lima.m.m2 <- inla(n ~ 1 + f(week,model = "rw1"),
                  control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
                  control.predictor = list( link = 1),
                  verbose         = F,
                  family          = "nbinomial",
                  data            =  db.n.masculino
                  )

lima.m.m3 <- inla(n ~ 1 + f(annus,model = "ar1") + f(week,model = "rw1"),
                  control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
                  control.predictor = list( link = 1),
                  verbose         = F,
                  family          = "nbinomial",
                  data            =  db.n.masculino
                  )
fit.m.m1 <- lima.m.m1$summary.fitted.values$mean
fit.m.m2 <- lima.m.m2$summary.fitted.values$mean
fit.m.m3 <- lima.m.m3$summary.fitted.values$mean
n.m      <- db.n.masculino$n


datos  <-  list("AR(1)"=fit.m.m1,"RW(1)"=fit.m.m2,"AR(1)RW(1)"=fit.m.m3) %>% 
           as.data.frame() %>% 
           melt( variable.name = "modelo",
                value.name    = "fit") %>% 
           group_by(modelo) %>% 
           mutate(actual    = n.m,
                  date      = db.n.masculino$date) %>% 
           ungroup() %>% 
           mutate(min_inter = c(lima.m.m1$summary.fitted.values$`0.025quant`,
                                lima.m.m2$summary.fitted.values$`0.025quant`,
                                lima.m.m3$summary.fitted.values$`0.025quant`),
                                
                                
                  max_inter = c(lima.m.m1$summary.fitted.values$`0.975quant`,
                                lima.m.m2$summary.fitted.values$`0.975quant`,
                                lima.m.m3$summary.fitted.values$`0.975quant`)  
           )


datos %>% 
ggplot()                        +
geom_line(aes(x=date,y=actual)) +
geom_line(aes(x=date,y=fit),color="#011f4b")    +
geom_ribbon(aes(x=date,y=fit,ymin=min_inter, ymax=max_inter), 
                alpha=0.4,       #transparency
                fill="#011f4b") +
facet_wrap(vars(modelo))        +
theme_bw()+
xlab(" ") +
ylab("Numero de casos")+
geom_vline(xintercept = as.Date("2019-01-07"), linetype="dotted", 
                color = "black", size=0.7)+
theme(plot.title = element_text(hjust = 0.5,size = "20pts"),
      strip.background =element_rect(fill="#011f4b"),
      strip.text = element_text(color = "white",face = "bold",size = "5pts")) 

2.2.1.2 Mujeres

#inla.list.models("likelihood")

lima.f.m1 <- inla(n ~ 1 + f(annus,model = "ar1"),
                  control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
                  control.predictor = list( link = 1),
                  verbose         = F,
                  family          = "nbinomial",
                  data            =  db.n.femenino
                  )


lima.f.m2 <- inla(n ~ 1 + f(week,model = "rw1"),
                  control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
                  control.predictor = list(link = 1),
                  verbose         = F,
                  family          = "nbinomial",
                  data            =  db.n.femenino
                  )

lima.f.m3 <- inla(n ~ 1 + f(annus,model = "ar1") + f(week,model = "rw1"),
                  control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
                  control.predictor = list( link = 1),
                  verbose         = F,
                  family          = "nbinomial",
                  data            =  db.n.femenino
                  )
fit.f.m1 <- lima.f.m1$summary.fitted.values$mean
fit.f.m2 <- lima.f.m2$summary.fitted.values$mean
fit.f.m3 <- lima.f.m3$summary.fitted.values$mean
n.f      <- db.n.femenino$n


datos  <-  list("AR(1)"=fit.f.m1,"RW(1)"=fit.f.m2,"AR(1)RW(1)"=fit.f.m3) %>% 
           as.data.frame() %>% 
           melt( variable.name = "modelo",
                value.name    = "fit") %>% 
           group_by(modelo) %>% 
           mutate(actual    = n.m,
                  date      = db.n.femenino$date) %>% 
           ungroup() %>% 
           mutate(min_inter = c(lima.f.m1$summary.fitted.values$`0.025quant`,
                                lima.f.m2$summary.fitted.values$`0.025quant`,
                                lima.f.m3$summary.fitted.values$`0.025quant`),
                                
                                
                  max_inter = c(lima.f.m1$summary.fitted.values$`0.975quant`,
                                lima.f.m2$summary.fitted.values$`0.975quant`,
                                lima.f.m3$summary.fitted.values$`0.975quant`)  
           ) 


datos %>% 
ggplot()                        +
geom_line(aes(x=date,y=actual)) +
geom_line(aes(x=date,y=fit),color="#011f4b")    +
geom_ribbon(aes(x=date,y=fit,ymin=min_inter, ymax=max_inter), 
                alpha=0.4,       #transparency
                fill="#011f4b") +
facet_wrap(vars(modelo))        +
theme_bw()+
xlab(" ") +
ylab("Numero de casos")+
geom_vline(xintercept = as.Date("2019-01-07"), linetype="dotted", 
                color = "black", size=0.7)+
theme(plot.title = element_text(hjust = 0.5,size = "20pts"),
      strip.background =element_rect(fill="#011f4b"),
      strip.text = element_text(color = "white",face = "bold",size = "5pts")) 

2.2.2 Performance

2.2.2.1 Hombres

perform.metrics <- metric_set(mae, rmse,msd)

tbl.yrd.m <- datos                    %>% 
             group_by(modelo)         %>%
             perform.metrics(truth    = n.m, 
                           estimate = fit)


tbl.yrd.m                               %>% 
pivot_wider(id_cols     = modelo,
            names_from  = .metric,
            values_from = .estimate)    %>%         
gt()                                    %>%
tab_header(
    title = md("Métricas de ajuste"),
    subtitle = "Hombres") %>% 
data_color(
    columns = vars(rmse,mae,msd),
    colors = scales::col_numeric(
      palette = c(
        "#011f4b","white"),
      domain = NULL)
  )
Métricas de ajuste
Hombres
modelo mae rmse msd
AR.1. 40.66904 50.09989 27.70859
RW.1. 34.64344 42.68089 27.66907
AR.1.RW.1. 30.38785 37.24076 27.67888

2.2.2.2 Mujeres

perform.metrics <- metric_set(mae, rmse,msd)

tbl.yrd.f <- datos                    %>% 
             group_by(modelo)         %>%
             perform.metrics(truth    = n.f, 
                             estimate = fit)



tbl.yrd.f                               %>% 
pivot_wider(id_cols     = modelo,
            names_from  = .metric,
            values_from = .estimate)    %>%         
gt()                                    %>%
tab_header(
    title = md("Métricas de ajuste"),
    subtitle = "Mujeres") %>% 
data_color(
    columns = vars(rmse,mae,msd),
    colors = scales::col_numeric(
      palette = c(
        "#011f4b","white"),
      domain = NULL)
  )
Métricas de ajuste
Mujeres
modelo mae rmse msd
AR.1. 30.25786 35.90613 0.10282570
RW.1. 24.32470 30.10790 0.06330288
AR.1.RW.1. 18.04854 22.13154 0.07311218

2.3 Prediccion: Enfoque Espacio-Temporal

El enfoque espacio-temporal, implica complementar el enfoque anterior complejizando el modelo tomando en cuenta la dimesion espacial. Para ello requerimos incluir efectos aleatorios espaciales para controlar por tal dimension de analisis. Cabe mencionar que modelos espaciales areales empleando INLA, requieren necesariamente la creacion de un grafo(graph) o una estructura de datos espaciales que captura la disposicion de las areas como una matriz de pesos espaciales.

\[\begin{align*} y_{t} = S_{t}+T_{t}+\nu_{t}+u_{t}+e_{t} \end{align*}\]

donde:

\(\nu_{t}\) : efectos estructurados

\(u_{t}\) : efectos no estructurados

Los 2 nuevos componentes \(\nu_{t},u_{t}\) pretenden tomar en consideracion en la prediccion la correlacion espacial presente en los datos . Al respecto :

  • Los efectos no estructurados corresponden a efectos aleatorios que pretenden controlar caracteristicas no observadas de cada area bajo estudio. Para modelarlas podemos emplear f(area, model = "iid")

  • Los efectos estructurados son efectos aleatorios que toman en cuenta explicitamente la estructura espacial

    • Pueden ser modelados en INLA de diversas maneras:
      • Empleando efectos espaciales besag :f(area, model = "besag")
      • Empleando efectos espaciales propios de besag :f(area, model = "besagproper")

2.3.1 Preprocesamiento

2.3.2 Ajuste y prediccion

2.3.2.1 Hombres

peru.m.m5   <- inla(n ~ 1 + f(annus,model = "linear") + 
                           f(week,model = "rw1")  + 
                           f(id.sp, model = "besag", graph = W.peru), 
                   data              = db.m.sp,
                   control.compute   = list(dic = TRUE, 
                                            waic = TRUE, 
                                            cpo = TRUE),
                   control.predictor = list( link = 1)
                  )

2.3.2.2 Mujeres

peru.f.m5   <- inla(n ~ 1 + f(annus,model = "ar1") + 
                           f(week,model = "rw1")  + 
                           f(id.sp, model = "besag", graph = W.peru), 
                   data              = db.f.sp,
                   control.compute   = list(dic = TRUE, 
                                            waic = TRUE, 
                                            cpo = TRUE),
                   control.predictor = list( link = 1)
                  )

2.3.3 Performance

2.3.3.1 Hombres

perform.metrics <- metric_set(mae, rmse,msd)

tbl.yrd.m <- datos                    %>% 
             group_by(modelo)         %>%
             perform.metrics(truth    = n.m, 
                             estimate = fit)


tbl.yrd.m                               %>% 
pivot_wider(id_cols     = modelo,
            names_from  = .metric,
            values_from = .estimate)    %>%         
gt()                                    %>%
tab_header(
    title = md("Métricas de ajuste"),
    subtitle = "Hombres") %>% 
data_color(
    columns = vars(rmse,mae,msd),
    colors = scales::col_numeric(
      palette = c(
        "#011f4b","white"),
      domain = NULL)
  )
Métricas de ajuste
Hombres
modelo mae rmse msd
AR.1. 30.25786 35.90613 0.10282570
RW.1. 24.32470 30.10790 0.06330288
AR.1.RW.1. 18.04854 22.13154 0.07311218

2.3.3.2 Mujeres

perform.metrics <- metric_set(mae, rmse,msd)

tbl.yrd.m <- datos                    %>% 
             group_by(modelo)         %>%
             perform.metrics(truth    = n.m, 
                             estimate = fit)


tbl.yrd.m                               %>% 
pivot_wider(id_cols     = modelo,
            names_from  = .metric,
            values_from = .estimate)    %>%         
gt()                                    %>%
tab_header(
    title = md("Métricas de ajuste"),
    subtitle = "Mujeres") %>% 
data_color(
    columns = vars(rmse,mae,msd),
    colors = scales::col_numeric(
      palette = c(
        "#011f4b","white"),
      domain = NULL)
  )
Métricas de ajuste
Mujeres
modelo mae rmse msd
AR.1. 30.25786 35.90613 0.10282570
RW.1. 24.32470 30.10790 0.06330288
AR.1.RW.1. 18.04854 22.13154 0.07311218